Light cone dynamics and reverse Kibble-Zurek mechanism in two-dimensional 

superfluids following a quantum quench 
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We study the dynamics of the relative phase of a bilayer of two-dimensional superfluids after the 
two superfluids have been decoupled. We find that on short time scales the relative phase shows 
"light cone" like dynamics and creates a metastable superfluid state, which can be supercritical. We 
also demonstrate similar light cone dynamics for the transverse field Ising model. On longer time 
scales the supercritical state relaxes to a disordered state due to dynamical vortex unbinding. This 
scenario of dynamically suppressed vortex proliferation constitutes a reverse-Kibble-Zurek effect. We 
study this effect both numerically using truncated Wigner approximation and analytically within 
a newly suggested time dependent renormalization group approach (RG). In particular, within RG 
we show that there are two possible fixed points for the real time evolution corresponding to the 
superfluid and normal steady states. So depending on the initial conditions and the microscopic 
parameters of the Hamiltonian the system undergoes a non-equilibrium phase transition of the 
Kosterlitz-Thouless type. The time scales for the vortex unbinding near the critical point are 
exponentially divergent, similar to the equilibrium case. 

PACS numbers: 



I. INTRODUCTION 

The technological advances of trapping and manip- 
ulating ultra-cold atom systems provide an opportu- 
nity to study many-body dynamics with unprecedented 
clarity. The realization of Bose-Einstein condensates 
in ultra-cold atom systems^, the Mott insulator transi- 
tion^, the BEC-BCS transition^, the Kosterlitz-Thouless 
transition^, demonstrated that this technology can be 
used as a quantum simulator of many-body phases. Here, 
the static state of a system in equilibrium is created 
and studied. Various dynamical aspects of ultra-cold 
atom systems have also been probed, such as dipole 
oscillations 7 , vortex excitations^, and soliton dynamics 9 , 
absence of equilibration in one-dimensional bosonic sys- 
tems^, spontaneous formation of vortices in spinor con- 
densatesii and many others (see Ref. [I4] for a recent 
review). In Ref. vortices excitations were created 
via laser stirring. In these experiments, the dynamics 
of only a few degrees of freedom were studied, such as 
the center of mass motion, or the dynamical evolution 
of a vortex. These experimental developments stimu- 
lated a considerable theoretical interest in understanding 
non-equilibrium quantum dynamics including analysis of 
dynamics following sudden quenches^, studying connec- 
tions between dynamics and thermodynamics 1 ^, dynam- 
ics through quantum critical points 16 . 

The focus of this paper is a detailed analysis of the 
full many body dynamics following the quench in a two- 
dimensional quantum rotor model. Physically we imag- 
ine the situation where two initially strongly coupled su- 
perfluids arc suddenly separated and we arc interested 
in the evolution of the relative phase between the two 
superfluids. In particular, we will be interested in the 
question of how the system relaxes to the equilibrium 



state. We note that experiments in a similar setup in- 
volving separation of two ID superfluids were reported 
in Ref. [2 21 ] and the corresponding theoretical analysis 
was done in Refs. [2314251] . Unlike the 2D case, phonon 
fluctuations in ID result in the exponential decay of the 
correlation functions and nonlinear effects in the form of 
phase slips do not bring qualitative changes to the be- 
havior of the correlation functions at least at low initial 
temperatures^. 

In equilibrium for the uncoupled layers, there are two 
possible phases. At low temperatures atoms in each layer 
(which we regard as identical) form a (quasi-)superfluid 
phase while at high temperatures they form a normal 
Bose gas. These phases can be distinguished by the long- 
range behavior of the single particle correlation function 
G(x) = (&t(o)6(x)> « p(exp[i(<£(x) - 0(0))]), where 6(x) 
is the single particle operator, p is the atom density, 
and 4> is the phase. We no te that a rotor representa- 
tion of bosons 6(x) ~ yj p(x) exp[i0(x)] is possible when 
the healing length characterizing the characteristic length 
scale of density fluctuations is short compared to other 
length scales in the problem. Under the same conditions 
the density fluctuations are negligible if we are interested 
in long distance physics. In the superfluid phase this 
function shows algebraic scaling, G(x) ~ |x|~ T / 4 , where 
the scaling exponent r is proportional to the temperature 
t ss T/T c , with T c being the Kosterlitz-Thouless temper- 
ature. At the transition point we have G(x) ~ |x|~ 1//4 . 
Above the transition, the correlation function shows ex- 
ponential scaling G(x) ~ exp(— |x|/£), with some corre- 
lation length £ which diverges near the transition tem- 
perature. The algebraic scaling of the superfluid phase is 
due the thermally excited phonon (Bogoliubov's) modes. 
In two dimensions these fluctuations generate quasi-long 
range order, rather than true long range order. The tran- 
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sition to the exponential regime is due to vortex excita- 
tions. Above the transition, vortex-antivortex pairs are 
deconfined so that vortices and anti-vortices become un- 
bound. These excitations generate a much more disor- 
dered phase field, which leads to exponential scaling of 
the correlation function. 

If we couple the two superfluids with a hopping term 
in the temperature regime of the critical temperature, 
the system forms a phase-locked state, see Refs. [20II21I ] . 
Here the correlation function of the relative phase scales 
as G(x) ~ exp(— |x|/£j) + C, where C^O, and is the 
correlation length of the phase-locked state. In this state 
the relative phase is well-aligned over long distances; its 
fluctuations are strongly suppressed. We then turn off 
the hopping and study the evolution of the system. 

In this paper we show that the relaxational dynam- 
ics occurs in two stages. The first fast stage, which we 
will term light cone relaxation, establishes a metastable 
quasi-equilibrium state of phonons (or Bogoliubov exci- 
tations) characterized by effective non-equilibrium tem- 
perature (in principle this metastable state can be com- 
pletely non-thermal). During this stage the correla- 
tions between two arbitrary points in space Xi and X2, 
G(xi,X2,£), where t is the time after the quench, ini- 
tially decay in time, independent of their spatial separa- 
tion x = |xi — x 2 | because these points are not causally 
connected: G(xi,X2,t) ~ l/t a , where a is a power-law 
exponent related to the parameters of the system. At a 
later time t*, when the condition 2vt* — x is fulfilled, 
these correlations (approximately) freeze in time so that 
G(xi,x 2 ,i) ~ \/x a . The exponent a thus defines the 
non-equilibrium phonon temperature in the system. Be- 
cause this first stage of dynamics involves only phonons, 
the exponent a can exceed the maximally allowed equi- 
librium value of one fourth, leading to a non-equilibrium 
super-critical metastable state, which can be thought of 
as a supercritical superfluid. It is analogous to an over- 
heated classical liquid, for which a liquid state can be 
sustained above the critical temperature if the creation 
of defects is avoided. We find that the power-law can be 
substantially above the critical scaling, and furthermore, 
that this metastable can be very long-lived. 

At longer time scales, vortex-antivortex pairs emerge 
and proliferate leading to the true equilibrium state. This 
process occurs at much longer time scales. We describe 
this thermalization process both numerically, using trun- 
cated Wigner approximation (TWA) and analytically. 
In particular, we show that thermalization (here cor- 
responding to the process of vortex-antivortex prolifer- 
ation) can be understood by extending renormalization 
group ideas to real time dynamics. By doing partial aver- 
aging over fast oscillating high energy degrees of freedom, 
we can rewrite the equations of motion of slower degrees 
of freedom through renormalized coupling constants. As 
in the case of equilibrium systems we observe two pos- 
sible scenarios corresponding to vortex-antivortex pairs 
being irrelevant (superfluid phase) or relevant (normal 
phase). Thus we are able to see how the system relaxes 



to one of the phases in real time. Divergent time (and 
length) scales in equilibrium systems translate into diver- 
gent relaxation times required to reach thermalization in 
the non-equilibrium case. 

Physically this decay of the metastable superfluid state 
to the new equilibrium is very reminiscent of the Kibble- 
Zurek (KZ) effect. The latter describes a ramp across 
a phase transition, starting on the disordered side. If 
the ordered state supports topological excitations, like 
vortices, then one expects very slow relaxation of the re- 
sulting state to the equilibrium due to vortex-antivortex 
recombination. This scenario is illustrated in Fig. Q] 
a): In the disordered phase we have excitations such as 
phonons, as well as topological defects. When we ap- 
ply a fast ramp across the phase transition, the phonon 
excitations thermalize on very short time scales, while 
topological defects can exist on much longer time scales. 
The mechanism of relaxation in our case is exactly com- 
plimentary and can be termed as reverse Kibble-Zurek 
effect. Here, the ramp across a phase transition starts 
from the ordered side, as illustrated in Fig. [T]b). In the 
ordered phase both phonon excitations and vortices are 
suppressed. When the system is ramped across the tran- 
sition, phonons are generated on a fast time scale. How- 
ever, vortices are generated at much longer time scales 
leading to the long-lived supercritical superfluid state. 
We point out that in thermally isolated systems (like cold 
atom systems) it is much easier to observe reverse KZ ef- 
fect because the disordered phase usually corresponds to 
a higher temperature. In isolated systems it is relatively 
easy to increase temperature by quenching some param- 
eter, while decreasing temperature requires much more 
effort and can be done only in open systems. 

This paper is organized as follows: In Sect. |H]we in- 
troduce the numerical method that we use and find that 
at short time scales the system shows light cone dynam- 
ics. In Sect. IIIII we consider the linearized dynamics of 
the bilayer system. Within this approximation both light 
cone dynamics and the emerging superfluid state can be 
understood. In Sect. IIVI we study the light cone dynam- 
ics of a solvable model, the transverse Ising model. In 
Sect. |V]we study dynamical vortex unbinding both with 
truncated Wigner approximation, and with a renormal- 
ization group approach. We note that a short version of 
this paper with some of the results was published ear- 
lier—. Here we expand the earlier treatment, present ad- 
ditional results and derivations, and formulate the real 
time renormalization group approach which explains the 
numerical results. 



II. MICROSCOPIC MODEL AND THE 
TRUNCATED WIGNER APPROXIMATION 

In this section we present the model that we use in our 
numerical approach. We consider two two-dimensional 
(quasi-)condensates that are aligned in parallel to each 
other, that are coupled by a hopping term which is then 
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FIG. 1: Illustration of the Kibble-Zurek (KZ) mechanism, 
which describes ramping across a phase transition from the 
disordered phase, and its counterpart, the reverse-Kibble- 
Zurek (rKZ) effect. The latter describes ramping across a 
transition from the ordered side. Its defining feature is the 
dynamical suppression of vortex unbinding, which happens 
on a much longer time scale than the appearance of phononic 
excitations. We propose to study the rKZ in a bilayer of 2D 
superfluids of ultra-cold atoms, by decoupling the superfluids 
and measuring the dynamics of the relative phase. 
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turned off. This can be achieved by increasing the po- 
tential between the two condensates. The coarse-grained 
Hamiltonian describing the relative phase (pi of the two 
superfluids corresponds to an XY model, to which we 
add a hopping term to describe the phase-locking in the 
initial state: 

H = fi ( - Y] -cos(0, - <f>j) + VV 2 

<ij> i 

-y(0^cos(V2^)), (i) 

i 

where is an overall (Josephson) energy scale, n de- 
scribes the ratio of kinetic and potential energies. We 
can formally replace these parameters by Q, k/tt = 2Jn, 
ttQo/k — U (so that Qq = s/2JnU, k = tt 
and V(t) = 2J±(t)n/H,Q, which gives a representation of 
two coupled Bose-Hubbard systems in the quantum rotor 
limited. In this limit the Bose operators are replaced by 
the phase-density representation and the fluctuations of 
density are assumed to be small. In the Bose-Hubbard 
model J is the in-plane hopping amplitude, U is the on- 
site interaction energy, n is the filling number, i.e. the 
number of particles per site, and Jj_ is the inter-layer 
hopping amplitude Jj_ . This representation gives at best 
a qualitative idea of how the model parameters relate 
to the parameters in experiment, but gives a more in- 
tuitive picture. We note that one can think about the 
continuum limit as discrete, where the lattice constant 
is approximately given by the zero-temperature healing 



FIG. 2: We simulate the dynamics of the relative phase of 
two 2D superfluids by solving the equations of motion and by 
averaging over the Wigner distribution of the initial state. A 
single run is shown here, for V = 100, k — 10 and T = 2, 
at the times t = 0, 5, 10, 20, 40, 100. Vortices are marked red, 
anti- vortices blue. 

length in the system, i.e. the length over which density 
fluctuations are suppressed. 

We emphasize that despite the BKT transition being 
classical in origin, i.e. driven by thermal fluctuations, the 
mechanism of vortex or phonon creation in the process 
we consider comes from quantum fluctuations. Indeed 
when the superfluids are strongly coupled together the 
density (which plays the role of momentum conjugate of 
the phase) strongly fluctuates because of the zero point 
motion. The heating mechanism of this system can be 
thought of as enhancement of this zero point motion fol- 
lowing the quench. 

It is convenient to introduce the rescaled quantities 
t = ftot/H, </> = \f^4>i an d n — \f\n. In terms of these, 
the classical equations of motion (EOMs) are 



dn, V 2, sr^ . (P(<pj ( - <pi)\ T .,,s„ . 7 /on, 

= -^-2^ sm (, ~^ j +V{t)psm0<pi{3) 

where we defined (3 — \/2tt/k. The indices 2% describe 
the four neighboring sites of site i. 
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We model the relative phase using a numerical im- 
plementation of the truncated Wigner approximation 
(TWA) (see Ref. [29| for a review): The expectation 
of any quantity at some time t > can be determined 
by sampling over a Wigner distribution at time t = 0, 
and solving the classical equations of motion from 
to t. This approximation is guaranteed to be accurate 
at short time o 30 ' 31 . This approximation is also exact 
for any quadratic theory so we expect it to be accu- 
rate in the first (light-cone) stage of dynamics primar- 
ily driven by phonon excitations. In our case we expect 
that TWA is also valid at longer times because when 
vortex-antivortex pairs start to emerge the system al- 
ready reached metastable state corresponding to finite 
effective temperature. At this point quantum fluctua- 
tions become suppressed by much stronger thermal fluc- 
tuations driving the slow vortex dynamics. 

We solve these EOMs for initial conditions that are 
distributed according to the Wigner distribution at t = 0. 
We can calculate this distribution under the assumption 
that J_l is larger than the other energy scales at t = 0. 
In this limit the phase fluctuations are small and can be 
described within the Bogoliubov approximation, where 
the system reduces to a sum of oscillators. The Fourier 
modes <f> q and h q at t = are distributed according to 
(see Ref. [H) 

/ <pi 2cr 2 n 2 \ 

where a = 1/ y/2u> qi r q ~ coth(w g /2To), and io q — 
v /4sin(g a; /2) 2 + 4sin(g y /2) 2 + V(3 2 with T being the 
initial temperature. Note that formally uj q diverges at 
V — > co. This divergence is unphysical, being an artifact 
of using Hamiltonian (fTJ) in the number phase represen- 
tation. In reality when J± becomes very large the trans- 
verse Josephson frequency saturates at w w 2Jj_. This 
happens at V ~ n or equivalcntly J± ~ Un. So for very 
strong initial coupling one can still use distribution (T3|) 
with V — > n. 

To visualize our simulations we show an example for a 
single run of the system on a 20-by-20 lattice in Fig. [5] 
The direction of the arrows on each lattice point describe 
the phase fa. We show 'snapshots' at various times. The 
plaquettes around which there is a phase winding of ±2tt 
are marked as vortices and anti-vortices. We see that at 
t = 0, the phases are well aligned due to the coupling be- 
tween the layers, with some small quantum fluctuations 
described by the Wigner function. The coupling is then 
turned off, vortices and anti-vortices are created pair- 
wise, and unbind on a long time scale, as we will discuss 
further on. To extract expectation values of our observ- 
ables from our simulations, we have to average them over 
many realizations of initial fluctuations. 

We use this method to extract the equal time correla- 
tion function: 

G(x,t) = (exp[»>/2&(t) - iy/24> j+x (t)]), (5) 




FIG. 3: Plot of short-time behavior of the correlation function 
as a function of time and space, at temperature T = 3, for k = 
10 and V = 20. The dynamics separates into instantaneous, 
damped oscillatory behavior, and a 'light cone' like pulse. 



where x is an integer separation between the points and t 
is the time after decoupling (see Fig. [3]) . Because we are 
using periodic boundary conditions G{x, t) depends only 
on the separation between the points x and does not de- 
pend on j. Note that this correlation function (or rather 
Jq dx'Gix' ,t)) can be directly measured in interference 
experiments ^ 19 i 22 . We indeed see very clear emergence of 
the light cone thermalization: At separations larger than 
2vt, where v is characteristic phonon velocity, G(x, t) is 
almost x independent - it uniformly decreases in time. 
Once 2vt > x the correlations freeze in time and de- 
pend only on x. The quantities in the system have been 
rescaled such that the phonon velocity is set to 1. 

We find that the state that emerges within the light 
cone shows algebraic scaling, and therefore can be re- 
ferred to as a superfluid. 



III. LINEARIZED DYNAMICS 

In this section we study the linearized dynamics of the 
system. Within this description, both the light cone dy- 
namics and the metastable SF state can be understood. 
The quadratic Hamiltonian describing the relative phase 
of two coupled superfluids reads: 

Ho = /^(-^(v«* + ^ + ^4 <*> 

v is the phonon velocity of the SF, approximately given 
by v — gn/ra. r$ is the short-range cut-off of the 
system, of the order of the healing length. J = v/r n 
is the KT energy. The term g±cj) 2 /2 is created by the 
hopping term of the bilayer. It is approximately given 
by g±_ — 4Jx?i. We note that when this Hamiltonian is 
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put on a lattice with a lattice constant n we obtain 



Hi — £7o 



(--£ 
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a , r o 2 
t>i + — n 



2r ; 



__(7) 

This expression can be also obtained directly linearizing 
the original Hamiltonian ([I]) . The index ji here describes 
the four neighboring sites of the site i, and is the 
filling fraction, related to the density via rij = nrf . Qq is 
related to the phonon velocity as f2o = v /r\ . We therefore 
find that the squeezing parameter k/tt is given by k/tt = 
r;/ro, i.e. it is the ratio of the discretization length scale 
ri and the short-range cut-off ro of the system, gx is 
related to V(t) by g x r? /2 = Q V(t). 

We now consider the time evolution of <p and n under 
([5]). It is convenient to go to momentum representation 
where different modes decouple from each other. Assum- 
ing also that we are interested in momenta smaller than 
1 jri , where the lattice effects are not important we obtain 
the following equations of motion: 



dt 
d 
dt 



-fi 



•( 



^0 



2V)6- k 



Sin — 71 



n 



(8) 
(9) 



where e? 



4 sin 2 k x /2 + 4 sin 2 k v /2, 



and k is dimen- 
sionless, k — —it... it. We rescale the time variable as 
t = Clot/H. The initial dispersion is then given by 



e 2 k + 2Vr /r l . 



J k.O 



(10) 



The dispersion u>). after the quench is simply u\ = e 2 , 
in these units. We solve these equations and calculate 
the equal-time correlation function at time t after the 
quench. We use 

G(x,t) = (exp(#(0,t))exp(-#(x,t))) (11) 
= cxp(-(<50 2 )/2), (12) 

where 5<f> = <p(0,t) — (f)(x,t). The averaging is now triv- 
ially done using the Wigner distribution (QJ. If we put 
the system back to the lattice we then find 



k 

' r k ,o 
^2u>k o 



(2-2 cos kx) 



cos 2 (wfci) 



sin 2 (w fe t) . (13) 



The quantities r^o and uiufi are defined as before. 

We now calculate the Green's function in the lin- 
earized regime numerically using Eqs. (fT2"j) and (jT3J) . We 
choose the discretization ri = ro, the initial tempera- 
ture T/Qq = 1) and the initial coupling Vf3 2 = 20. In 
Fig. @] we plot (8(f) 2 ) and in Fig. [5] we plot the correla- 
tion function. In both plots the light-cone dynamics is 
clearly visible. Because of the translational invariance 
the correlation function that emerges in the light-cone 
only depends on the relative distance is given by 



G(x,t)«C 1 |x|- T */ 4T ^ 



(14) 




200 



2 ) of the linearized system, for T = 1 and V/3 2 



FIG. 4: 

20, as function of the lattice site, and vt 
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FIG. 5: The correlation function of the linearized system, for 
T = 1 and V/3 2 — 20, as function of the lattice site, and vt. 



for x <C 2vt, where T* is an effective temperature that is 
estimated below, and C\ is a numerical prefactor. Out- 
side of the light cone (x 3> 2vt) the function G(x, t) only 
depends on time t but not on the distance x: 



G(x,t) = C 2 \t\ 



-T*/4T K 



(15) 



where T* is the same effective temperature. At the light 
cone boundary x ~ 2vt the two asymptotics for the cor- 
relation function l|14|) and (|15p approximately coincide. 
However, we note that the prefactor G2 is in general dif- 
ferent from C\v~ t *I ATkt as it is evident from the exis- 
tence of a wavefront that is visible in Figs. 2]and[5j 

The temperature that emerges inside the light cone can 
be estimated by considering the quadratures of 4> at long 
times: 



rk.o 



4ui 



k.O 



r k ,oUk,o 
*4 ' 



(16) 



We find that the whole Wigner function in the non- 
interacting evolution remains Gaussian. It means that 
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for each mode the Wigner function is equivalent to that 
of a harmonic oscillator at finite 'temperature' T£ , which 
is in general mode-dependent: 



u k 



2ul 



where r£ = 1/ 

n = 



2uj k fl 

= l/tanh(cj fe /2T*). Solving for T* gi 



;ives 



2 tanh 



-l 



'^^tanh( WM /2T) 



(17) 



(18) 



For large V/3 2 , which corresponds to initially strong cou- 
pling between two superfluids, this simplifies to a single 
value, independent of k: 



4tanh(VW?2/2T)' 



(19) 



For small initial temperatures T we have T* fa 0^/4 
(T* = 2J±/J in terms of the original Hubbard param- 
eters), that is, the temperature is fully determined by 
the initial coupling energy. The coupling energy be- 
tween the two layers is transferred into the in-plane ki- 
netic energy. We remind again that this result is valid 
as long as J± < Un, otherwise the dependence of T* on 
J_l saturates and for the infinite coupling limit we have 
T* ~ Un/ J. For large T we have T* pa T/2. This result 
is a reflection of the doubling of the degrees of freedom 
when two layers are uncoupled. In Fig. [5] a) we show 
dependence T*(T) evaluated according to Eq. (fl"9jl . for 
V = 20, and for k = 1,3,10 corresponding to lowering 
J±. For k = 1, T* is always above the critical tempera- 
ture T c = 7r/2, for k = 10, it crosses it. We therefore ex- 
pect to see very little vortex formation for small temper- 
atures for k = 10, and many vortices for all temperatures 
for k = 1. The intermediate value k — 3 approximately 
describes the transition between these limits. To make 
this point more clear in Fig. [5] b)-d) we plot full non- 
linear TWA simulation for each of the three cases. We 
run the quench for different temperatures T, and plot 
the number of vortices n v in the system as a function 
of time. This number is obtained by counting the vor- 
tices (indicated by red plaquettes in Fig. [2J, and then 
by averaging over many runs. We find that for k = 1 
the number of vortices is virtually unchanged implying 
that the dynamics is completely dominated by quantum 
fluctuations, whereas for k = 10 this number drops to 
zero when the temperature is lowered. These results are 
consistent with the emergent temperature T* obtained 
within the linearized approach. 



IV. LIGHT-CONE DYNAMICS IN THE 
TRANSVERSE ISING MODEL 

In this section we demonstrate that light cone dynam- 
ics is not just characteristic for the system we are inter- 
ested in, which is characterized by low energy bosonic 
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FIG. 6: a) T* , as given in Eq. Q2] for n = 1,3, 10, from top 
to bottom, and for V = 20. The line T c = n/2 was added 
to indicate the critical temperature, b) - d) Simulations for 
these values of k and V. We plot the number of vortices n v 
as a function of time t and initial temperature To. 



wave excitations. The same mechanism of reaching 
a steady state is much more general and is likely re- 
lated to the existence of the maximum group velocity in 
Schrodinger systems as was proven by Lieb and Robin- 
son^. In this section we demonstrate the presence of 
the light-cone dynamics in another solvable model, the 
transverse Ising chai n 26 ' 33 described by the Hamiltonian 



(20) 



where Jj is an overall energy scale, g describes the 
strength of the transverse field, and a x,z are the Pauli 
matrices. We follow the calculational procedure in 
Ref. [33j]. First, we use a Jordan- Wigner transformation 



1 - 2n, 



11(1-2^)^ + 4), 



(21) 
(22) 



where Ci are Fermi operators, and = c\ci. This 
transformation leads to a fcrmionic representation of the 
Hamiltonian, that can be further diagonalized using the 
Bogoliubov transformation 



U ky gC k - iVk,gC*_ k 



(23) 



where c k is the Fourier transform of Cj, u k g and v kl9 
are given by cos(#fc iff /2) and sm(9 kyg /2), where 9 kt9 = 
arctan(sin k/(g — cosfc)). The resulting dispersion is 



2,7/vV - 2# cosfc + 1. 



(24) 
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FIG. 7: The correlation function (eg (t)cr^(i)) for a quench 
from g = 3 to g' = 1, in (a), and from g = 3 to </ = 0.5, in 
(b), as a function of time t, specifically of 2 Jt, and the spatial 
distance r. 



We consider a time dependent g(t). For t < we have 
= g, and we assume the system to be in equilibrium. 
We then assume that for t > 0, g(t) jumps to the value 
g' . The equal-time correlation function of of can be cal- 
culated exactly by expressing it in terms of the operator 
th, i.e. (af(t)olj(t)} = 1 — 4(nj(t)) + k{ni(t)rij (t)). It can 
be shown that the average density fermionic density (cor- 
responding to the z-component of the magnetization) is 
given by 



(m(t)) = MO)) 



>{k,t) (25) 



with 



(cos(2e fc , g ,t/fc)/2 - 1/2) {g' - g) sin 2 k 



y/g 2 - 2gcosk + l(g' 2 - 2g' cos k + 1) 



(26) 



and 



K(o)> 



i 



i 



\ - 



g — cos k 



(27) 



2 2M ^ ^. g 2_2 5 cosfc + l' 
In turn the density-density correlation function 
(rii(t)nj(t)) reads 



= ^ E ( ex P(- l ( fc i- fc 2)(n-r,))((i;2 iig+ F^Kfci^))(< :S -^, S '(fc2,t)) 
+ K Il9 ^ ll9 + Gg,g'(ki,t))(uk 2 , g Vk 2 ,g + G; )fl ,(Afc,i))) + «, 3 +F ffiff ,(fci, *))«,„ + ^,s'(fc2,i))), (28) 



where 



Gg, g >(k,t) = (iBin(2e k , g t/h) + 



1 (cos(2ej. .gd/h) — l)(g' — cosfc) 

2 Vff' 2 - V cosfc + 1 



(.9' - 5) sin fc 



^{g 1 - 2gcosk + l)(g' 2 - 2g' cosfc + 1)< 



(29) 



Using these expressions we can easily analyze the quench 
dynamics. In Fig. [7] we show two examples showing spin- 
spin (density-density) correlation function after a quench. 
The first example corresponds to the ramp from g = 3 to 
g' = 1, i.e. a quench to the quantum critical point. At 
this point the dispersion (f24| becomes gapless and linear 



at small energies. Then the light cone dynamics is an- 
ticipated because there is a well defined "speed of light" 
characterizing the propagation of excitations, which is 
equal to 2 J. Indeed Fig. |7K) shows clear signature of 
such dynamics. There is a clearly visible "light cone", 
which separates into an instantaneous part (connecting 



not causally connected points) that is independent of the 
distance, and a spatially dependent (causal) part, that 
expands in the form of a wave front. In Fig. b) we use 
g' = 0.5, where at low energies the spectrum of excita- 
tions is gapped. Although the dispersion in this model 
is linear (relativistic) only at sufficiently high energies 
above the gap we still see a clear light-cone structure. 
The expansion velocity of the "light cone" in this case is 
consistent with twice the maximum of the group velocity 
v gr (k) = de/dk given by 

f 2Jg for \g\ < 1 

V gr ,rna X ~ ^ 2 J for \g\ > 1. [ ^ > 

So for g' — 1 we find 2v™. TOOa . = AJ, and for g' = 0.5 
we find 2v gr)max — 2J. These are indeed the expansion 
velocities that we see in Fig. [7] 



V. DYNAMICAL VORTEX UNBINDING 

In this section we address the important question of 
how the supercritical state relaxes to the ground state, 
i.e. the second stage of the dynamics. As we mentioned 
in the introduction the anticipated mechanism for this re- 
laxation is vortex unbinding. This process is intrinsically 
nonlinear and requires a more sophisticated treatment 
than that of the noninteracting "light cone" dynamics. 
In this work we use two complimentary approaches. In 
Sec. IV Al we use a numerical implementation of the TWA 
to simulate the dynamics in the system. In Sec. IV Bl we 
generalize a renormalization group approach to analyti- 
cally describe the process of relaxation in real time. 



A. Numerical approach 

Within TWA we need to solve the full nonlinear equa- 
tions of motion (j3J) subject to the initial conditions dis- 
tributed according to the Wigner function (j4| . Then the 
equal-time correlation functions or other observables are 
found by averaging the Weyl symbol of the corresponding 
observable computed at time t over the fluctuating ini- 
tial conditions. Note that since we are interested only in 
phase-phase correlation function the corresponding Weyl 
symbol is obtained by simply substituting the Heisenberg 
quantum operator corresponding to the phase with the 
classical phased. In Fig. |8] we show the result of such 
simulations. We can observe how the metastable super- 
fluid state relaxes to the disordered state. For that, we 
show the correlation functions of the system on a much 
longer time scale than in Fig. [3l The exponent of the al- 
gebraic scaling gradually decreases. Eventually the cor- 
relation function is more accurately approximated by an 
exponential fitting function, signalling that the thermal 
Bose gas phase has been reached. Because this is the 
phase of deconfined vortices, and because the intermedi- 
ate superfluid phase is well described by a phonon-only 




FIG. 8: Long-time behavior of the correlation function for 
T = 1, k = 8, and V = 80. The correlation function first 
develops algebraic scaling, so the system forms a metastable 
quasi-superfluid state. On longer time scales the correlation 
function shows exponential decay. The coherence is lost due 
to dynamical vortex unbinding. 

description, we conclude that the dynamical transition 
that we observe is due to vortex unbinding. The exam- 
ple of a single run shown in Fig. is consistent with 
this picture: Defects are created soon after the quench, 
but they only gradually separate on a much longer times 
scale. It is this process that we refer to as the reverse 
Kibble- Zurek mechanism. 

To better characterize the process of vortex unbinding 
further we fit the correlation function G(x,t) to either 
algebraic or exponential fitting functions. Such choice is 
motivated by the two possible regimes of the equilibrium 
system and is supported by the analytic renormalization 
group results presented in the next section. The alge- 
braic fitting function we use is c(L/7r| sin(7ra:/L)|)~ T / 4 
and the exponential function is ccxp(— | sin(7rx/L)|/a;o). 
Note that in the fitting functions we use the conformal 
distance L/n\ sin(7ra;/L)|), which is more appropriate in 
finite systems with periodic boundary conditions (see e.g. 
Ref. [33j|). In equilibrium the algebraic exponent r would 
be the relative temperature T/T c . Any value above 1 
is therefore supercritical. The parameter xq defines the 
length scale of the exponential decay. The parameter c 
in both functions gives an overall scale. 

Using these fitting functions we analyze four different 
situations corresponding to the same initial temperature 
T = 1 and the same parameter n = 8, but with dif- 
ferent initial couplings V between the planes. The first 
(I) case corresponding to V = 80 is identical to the one 
plotted in Fig. |8] The other three curves correspond 
to V = 70,50,20 (II-IV). In Fig. M we show the expo- 
nent t extracted from the fit as a function of time for 
these situations. In all of them at short times G(x,t) 
develops algebraic scaling when the light-cone dynamics 
reaches the system boundaries. For the cases I— III the 
emerging scaling exponent r is well above the critical ex- 



9 




FIG. 9: Time dependence of the exponent r extracted from 
fitting the long-time correlation function G(x,t), for different 
initial couplings. In all four examples we use T = 1 and k = 8. 
The initial couplings V are chosen as V = 80, 70, 50, and 20, 
corresponding to curves I to IV. The curve I corresponds to 
the example shown in Fig. [8] In this case the correlation 
function can be well fitted with an algebraic function for up 
to t ~ 60, after that G[x, t) is better fitted by an exponen- 
tial function with a decay length of the order of the lattice 
constant. For the other cases, G(x, t) is well fitted with an 
algebraic function throughout the whole time interval. 



ponent. After that, the exponent gradually increases on 
much longer time scales. During this process, the decay 
of the correlation function is still fitted well with the alge- 
braic function. Eventually the algebraic scalings breaks 
down and G(x, r) develops exponential scaling, indicat- 
ing vortex unbinding. This regime of exponential scaling 
is reached for V — 80 (I) within the time interval shown 
in Fig. H For V = 70 (II) and V = 50 (III) the time 
scale of the vortex unbinding is longer then the time in- 
terval shown. For V = 20 (IV) the system equilibrates to 
the superfluid state. Because in this case the exponent r 
is less than one, vortices never unbind and the algebraic 
scaling persists at all times. We conclude from these ex- 
amples that there can be a sizeable range of initial values 
of V which generates the scenario of a supercritical super- 
fluid, and of dynamically suppressed vortex unbinding. 
Furthermore, the algebraic scaling exponents that can 
occur in the metastable state are well above criticality, 
and should be easily distinguishable from subcritical val- 
ues. These supercritical exponents can be detected using 
interference experiments along the lines of Refs. I5IIT9I. 



B. Renormalization group approach 

In this section we develop the renormalization group 
(RG) approach to dynamical vortex unbinding. We find 
that the dynamical evolution of the system can be re- 
lated to the RG flow of the equilibrium system. The idea 
of RG in real time is quite similar in spirit to the RG in 
imaginary time. Namely our goal is to eliminate high en- 



ergy, high momentum degrees of freedom. In equilibrium, 
this is done by the means of usual perturbation theory 
(or Gaussian integration), which is justified because of 
the large energy gap separating high energy states from 
the low-energy degrees of freedom we are interested in. 
In real time the idea of renormalization is quite similar. 
High energy (momentum) phonons are not very sensi- 
tive to slow processes leading to vortex formations. Thus 
these phonons can be well treated within the linearized 
approach. However, due to nonlinearities such phonons 
slightly renormalize the parameters governing dynamics 
of low energy degrees of freedom. This renormalization 
is precisely what we are interested in. Note that techni- 
cally in the RG procedure we perform averaging of the 
equations of motion over short times. Then odd powers 
of highly oscillating fields average to zero while averag- 
ing of the even powers gives some constant contribution. 
This contribution is precisely what renormalizes coupling 
constants governing the low temperature dynamics. 

We point out that typical RG flow diagrams contain 
mostly non-equilibrium points, in fact, all except for 
the fixed points. As we have seen in the previous sec- 
tions, one can associate an effective temperature to the 
metastable state that emerges after the dephasing of the 
phonon modes. In turn with this effective parameter we 
can associate a location of the transient state in the RG 
flow of the equilibrium system. This effective tempera- 
ture can then either gradually increase, until the system 
starts to show exponential scaling, or the system can al- 
ways remain superfluid, if the algebraic scaling is subcrit- 
ical and the effective temperature always remains below 
Tkt- This behavior resembles the equilibrium RG flow 
of a Kosterlitz-Thouless transition (which now occurs in 
real, not imaginary, time) , on which we elaborate in this 
section. 

Instead of directly analyzing the rotor model to de- 
scribe the Kosterlitz-Thouless physics and vortex unbind- 
ing, we will work with the dual Z\ clock model (or equiv- 
alently 2D sine-Gordon model), described by the action 

S = J d 2 r(^{d x 9) 2 -^cose). (31) 

For the details of the duality transformation see Ref. [l7| • 
The parameters of this model can be related to those of 
the XY model by 

IT IT 32 
8ir Tkt 47t 2 Jkt 
| = expt-Sy, (33) 

where E c — S C T is the vortex core energy, and A is a 
measure of the relative temperature. We note that the 
action in Eq. [31] has a high-momentum cut-off A, which 
is the inverse of the short-range cut-off a, i.e. we set 
Aa = 1. To describe the dynamics of this model we use 
the effective 2D sine Gordon Hamiltonian: 

H/T = Jd\(^-±(d x 6f + ^cos6). (34) 
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Here the parameter fi is chosen, so that the dispersion of 
the linearized XY model is recovered: 



Xk 2 T 2 ' 



(35) 



which we also write as lou = v\k\, where the velocity v 
is given by v = \J /jAT 2 . The nonlinear term cos 9 in 
Eq. (f34|) describes the vortex field. If this term is im- 
portant (large g) then the field 9 localizes corresponding 
to a highly disordered phase of the dual field tp, i.e. to 
the normal state. Conversely small g corresponds to the 
superfluid algebraic regime. The starting point of our 
RG analysis will be supercritical superfluid state which 
emerges after short time light cone dynamics. Because 
the Kosterlitz-Thouless transition is classical in nature 
occuring at high temperatures the quantum fluctuations 
are no longer expected to be important and instead of 
Wigner function as the new initial condition we can use 
its classical Boltzmann's limit. The initial state for the 
vortex dynamics, described by the effective temperature 
T, is thus fully characterized by the quadratures of the 
spectrum: 



(Pkpk) 



1 

Xk 2 ' 
1 _ 



XT 



(36) 
(37) 



The equations of motion corresponding to the Hamilto- 
nian (|34[) are given by 



dt 1 
dt 



a 1 

fiTp. 



(38) 
(39) 



We now apply the following renormalization procedure 
to these equations. We rescale the spatial and temporal 
variables as r -» r(l + dA/A) and t — > t(l + dA/A), and 
the p-field as p — > p( 1 — dA /A) . This implies that the mo- 
mentum cut-off A is rescaled as A — > A' = A(l— dA/A), so 
the momentum degrees of freedom between A' and A are 
removed. Without the non-linear term in Eq. I3"51 these 
rescalings leave the equations of motion invariant. The 
linear dynamical evolution can therefore be considered to 
be the non-interacting fixed point of the RG. We now ask 
the question, how this dynamical evolution is affected by 
the non-linear term. Specifically we want to determine 
how the equations of motion behave at long times and 
distances. For this, we go beyond the bare rescaling and 
correct for the integrated-out degrees of freedom up to 
second order in g. The resulting flow equations are of 
the well-known BKT form: 



dg 
dl 
dX 

~dl 



2-— )g 



(40) 
(41) 



dl = dt/t = dr/r 




-a q H(r,t,A,g) 
G>pH(r,t,A,g) 



p(r',t') = -d q H(r',t',A',g') 
e(r',t') = S p H(r',t',A',g') 



FIG. 10: Schematic representation of a renormalization step 
in the real-time RG approach. In each step we renormalize 
simultaneously the space and time variables. This 'moves' the 
equations of motion from (r, t) to (r',t'). We correct for the 
integrated-out degrees of freedom to second order in g, which 
renormalizes the parameters g and A according to Eqs. 1401 
andSU 



where I = In A, and a is a non-universal prefactor. The 
RG step generated the equations of motion at time t' and 
distance r' from the equations at time t and distance r, 
with renormalized coefficients, according to Eqs. 0D] and 
W\\ Therefore the time dependence of the coefficients can 
be read off the solution of the RG flow, by realizing that: 
dt/dl = t or t = t Q e l . In Fig. [TU]we show a schematic 
representation of our RG process. In the Appendix we 
discuss the derivation of the flow equations and give their 
more complete form. 

One conclusion from Eqs. 1301 an d EH is that the crit- 
ical exponent of the dynamical process is equal to the 
one of the equilibrium system. We see from Eq. 00] that 
the critical value of A is A c = 1/87T, which corresponds to 
T = Txt as can be seen from Eq. 1321 Another important 
observation is that the RG equations (14"01) and (l4"Tj) pre- 
dict a non-equilibrium analogue of the BKT transition, 
where depending on the initial fluctuations in the sys- 
tem, the vortex-antivortex pairs can either unbind in the 
long time limit or remain bounded. This transition, as 
in the equilibrium case, is characterized by exponentially 
divergent time and length scales. Physically these diver- 
gencies correspond to a very slow process of equilibration 
of vortices near the nonequilibrium phase transition. 

We can also use the RG flow to determine the time 
scale of vortex unbinding by using 



g(t*) ~ 1. 



(42) 



When T is well above Tkt, the time scale can be deter- 
mined from Eq. 1401 



t* 



exp 



( E c /2 \ 
\T-T KT ) 



(43) 



where E c = S C T . Away from the transition, the time 
scale of vortex unbinding is therefore exponentially in- 
creased, because of the energy cost given by the vortex 
core energy. Very close to the transition t* scales as: 



t* ~ exp(exp(-S c /2)/v/l - T KT /T). 



(44) 
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The time scale is renormalized because of the critical scal- 
ing in the vicinity of the transition. 

VI. CONCLUSIONS 

In conclusion, we have studied the dynamics of the rel- 
ative phase of a bilayer of superfluids in 2D, after the hop- 
ping between them has been turned off rapidly. We find 
that on short time scales the dynamics of the correlation 
function shows a 'light-cone'-like behavior. Depending 
on the parameters of the system, the light cone dynam- 
ics can result in a phase that shows supercritical algebraic 
scaling, and can therefore be thought of as a superheated 
superfluid. On long time scales the system relaxes to a 
disordered state via vortex unbinding, which constitutes 
a reverse-Kibble-Zurek mechanism. The properties of the 
dynamical process can be understood with a renormal- 
ization group approach. We find that the dynamical evo- 
lution of the system resembles the RG flow of the equi- 
librium system. In particular, using the RG equations 
we found two possible scenarios of the system reaching 
the steady state: (i) if initial quantum and thermal fluc- 
tuations are weak the vortices are irrelevant and long 
time long distance behavior is governed by the algebraic 
fixed point. The only role of vortices is then renormal- 
ization of the superfluid stiffness and the sound velocity, 
(ii) If the initial fluctuations are strong then the vor- 
tices become relevant and proliferate resulting in a nor- 
mal (non-superfluid) steady state. In this case RG gives 
the time scale of vortex unbinding, which exponentially 
diverges as the system approaches the non-equilibrium 
phase transition. The behavior of the relative of phase of 
two superfluids can be accurately studied by interference 
experiments of ultra-cold atom systems, and therefore 
our predictions are of direct relevance to experiment. 

Acknowledgments 

We thank A. Castro Neto for useful discussions. Work 
of A. P. was supported by NSF DMR-0907039, AFOSR, 
and Sloan Foundation. L.M. acknowledges support from 
NRC/NIST, NSF Physics Frontier Grant PHY-0822671 
and Boston University visitor's program. 

Appendix A 

In this Appendix we derive the RG Eqs. 00] and EJ 
which can also be written as a second order differential 
equation for 9 

~4* 9 = AA0+4sin#. (Al) 
\x d,t z a 1 

To simplify the derivation, here and throughout the Ap- 
pendix, we formally change notations AT — > A, /iT — > /i, 



and gT — > g. The idea of momentum shell RG is that we 
treat high momentum components of 9 and p (or equiv- 
alently 9) perturbatively, while not making any approx- 
imations about the low momentum components. Our 
goal is to find renormalization of the equations of motion 
governing the low momentum components. So we split 

9(r,t) = e<(r,t)+9>(v,t) 7 (A2) 

where the Fourier expansion of 6 > (r, t) only contains mo- 
menta in the shell A' = A - SA < \k\ < A and 9 < {r,t) 
contains all other Fourier components: 

9<(r) = -= V exp(ikr)6» k (A3) 
9>{v) = -JL Y, exp(»kr)0 k . (A4) 

V V A'<fc<A 

We will treat 8 > (and correspondingly p > ) pertur- 
bartively in g since the nonlinear term should only weakly 
couple to the high frequency field. We expand the high- 
momentum field as 

9 > (k,t) = 9>(k,t) + 6>(k,t). (A5) 

Here 9q (k, t) is the solution of the equations of motion, 
with g set to zero: 

9>(k,t) = -%> k sin(w A i) + 0jf k cos(w A *), (A6) 

where uik = v\k\, and the velocity v is v — y/Xjl. In the 
next leading order we have 

0>(k,i) = J drF^k, t) sin(u> A (* - r)), (A7) 
where 

F x {k,r) = J d 2 rexpHkr]sin(6l<(r,T)), (A8) 

and we used Aa = 1. Note that in the last equation in 
the argument of the sinus we changed 9$ to 9q because 
the contribution from 9q is smaller by the factor SA/A. 
So we see that in the leading order in g the high mo- 
mentum component of 9 oscillates with time at very high 
frequency uja ■ In the next order in g the high momentum 
component also acquires a low frequency component (as 
we will discuss below). 

Next we consider the equation of motion (|A1|) expand- 
ing it up to the second order in 9 > : 

~^{v,t) « AA0(r,i) + Acos(0<(r,t))0>(r,i) 

Because of the nonlinearity high-momentum modes cou- 
ple to the low momentum modes leading to the renor- 
malization of the couplings governing the dynamics of 
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the latter. The idea of RG is to average equations of 
motion for low-momentum (slow) components over the 
fast oscillations. The averaging is trivially done in the 
last term of Eq. (|A9|) . There it is sufficient to use zeroth 

order in > . Using that sin 2 (wA£), cos 2 (oj\t) — 1/2 we 
find that averaging of the last term simply renormalizes 
the coupling g: 

where Ek is the average energy o f the mode k over the 
period (we used the fact that Xk 2 \8k\ 2 = Ek)- We note 
that in a Boltzmann ensemble, we would have E^ = 1, 
because the energies here are in units of the tempera- 
ture T. With this assumption we would recover the flow 
equation of the equilibrium case. 

Instead of this assumption, we proceed by noting 
that under RG transformations coupling constants slowly 
change in time. This implies that the adiabatic invariants 
per each mode are approximately conserved, as discussed 
in Ref. [34]. For an oscillator the adiabatic invariant is 
Ik = Ek/^k- Thus we see that the energy of the mode 
is proportional to the frequency. Noting that at initial 
time Ek(t = 0) = 1 (in non-rescaled units this would be 
Ef.(t — 0) = T , where T is the initial non-equilibrium 
temperature), one can rewrite Eq. (|A10p as follows: 

where we introduced the analogue of the Luttinger- 
Liquid parameter K = \J X/fi. v is the velocity which 
is now given by v — vA/U (note that in the original, not 
rescaled units, v = T^/Xfl). 

Next let us consider the second term in Eq. (|A9|) . This 
term is more subtle since if we use 0q the average over 
fast fluctuations will give zero. So we need to use the 
first correction 0> , which would be a correction at sec- 
ond order in g. We note that the linear term in Eq. (|A9|) 



can also be expanded to second in g, generating simi- 
lar contributions. However, when written as Eq. (|A1|) . 
such a term a cancelled by a corresponding term from 
expanding d 2 8/dt 2 . Alternatively we can view Eq. (IA9|) 
as written for the < component, this automatically en- 
sures that only the nonlinear term is responsible for the 
renormalization. 

Let us look closer into the Eq. (|A7[) . We are dealing 
with the integral over the fast oscillating function of r: 
sin(wA(£ — t)) and the slow oscillating function F. This 
integral can be evaluated by integrating by parts: 

/* **i(r) «dn( WA (t - r)) = F l( r) C0S ^ A( *- r)) * 

JO W A 

/ dr \ K ' cos(u} A (t-T)). (A12) 

wa Jo dr 

Note that the second integral contains a large denomina- 
tor 1/uja- In the first term only the limit t — t gives a 
nonoscillating contribution to the integral. We can con- 
tinue the expansion in powers of 1/u>a- Note that the 
next term proportional to F\ will contain only highly os- 
cillatory part and can be neglected. So up to the third 
order in 1/oja we find: 

f* , / x / , \ \ FAt) 1 d 2 F 1 (t) , K . 
/ drFx(T sin {uj A (t-T )« 3—^. A13 

Combining Eqs. ijAT]) . (|A8)l . and (lA13|) we find 

0({k,t) w I / d 2 rexp[-ikr] sm(0<(r,t)) 
--^2 / d 2 rexp[-ikr]cos(6» < (r,i))6'' < (r,t)(A14) 

where we used Aa = 1 again. Here we neglected by the 
term proportional to (9 (v, t)) 2 , because it leads to a 
subdominant (in the RG sense) contribution. From this 
we find that 



*r(r,*)«f / ^e^l^e-^sin^Cx,*))-^^ J ^e lkr |d 2 xe- kx cos(0<(x,i))0<(x,t). (A15) 

shell shell 

We now consider the term (g/a 2 ) cos(^ < (r, £))0^"(r, t): 

^cos(0<(r,t))0Ur,t) « £2 J (S2 eikr / d 2 xe- ik *±sm(8<(x,t)-9 < (r,t)) 

shell 

/ ^/'"■*"5 i< « I - t »' (A16) 

shell 

where we neglected terms such as sin(20 < (r, t)). Because we are integrating over high momentum shell this integral 
will be suppressed unless x is close to r. This suggests change of variables x = r + £ and Taylor expanding < (x, t) 
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in powers of , i.e. 0<(x, t) « 0<(r, t) + £V0<(r, t) + f&*£/j§J^. Then 



where 
A 2 



d 2 £Jo(A£), 



C2 
A 4 



| d 2 a 2 J (AO. (A18) 



When we need to substitute these expressions back into 
Eq. (|A9p . we find that the term containing AO renormal- 
izes the coupling A as 



A — !• A + 



C 2 g 2 SA 
fell' 



(A19) 



In addition there is an extra term proportional to 9 gen- 
erated in Eq. (|A9[) . which renormalizes /i: 



1 



1 Ci g 2 SA 

fl 47T A 2 /Lt A ' 



(A20) 



Finally we restore the cutoff by rescaling k —> k(l — 
SA/A), r -> r(l + <5A/A), t -> t(l + 5A/A), and p -> 
p(l — 6 A/ A). This rescaling additionally renormalizes 
the coupling g: g — > g(l + 25A/A). Combining this re- 
sult with Eqs. (jAlTj) . [|Al"9"j> . (|A^0|) we find the following 
renormalization group equations: 



dl 
dK 

~dT 

dv 



g(2 
1 .</ 



1 1 



16^7^ v 2 

r.2 



Attvo K / 
2 (C 2 + 2d), 
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16ttvK 2 



(C 2 - 2d). 



(A21) 
(A22) 
(A23) 



where Z = In A. We can read off from Eq. IA231 that 
if the system contains a fixed velocity, for example in 
relativistic systems, we need to have C2 = 2Ci, to enforce 
that the velocity is invariant under the flow. 

Note that if the initial system is already close to the 
critical point then the RG equations above simplify to 



<ig_ 

dl 
dX 
dl 



g 

8tt A ' 



1 

4^rA 



(A24) 
(A25) 



which are equivalent to Eqs. (|40| and (|4Tj) . Note that a 
more complete set of RG equations (|A21I) - (|A23|) has the 
same universal predictions of the dynamical phase tran- 
sitions and exponential divergence of the time scales as 
the simplified equations above. Also note that the real 
RG equations bear close analogy to the flow equations in 
imaginary time characterizing the equilibrium Kosterlitz- 
Thouless transition^. Thus the non-equilibrium KT 
transition discussed here is characterized by exponen- 
tially divergent length and time scales. Physically these 
long scales characterize very slow process of vortex un- 
binding and equilibration at long distances. Note that 
the RG equations (|A21[) - (IA23[) also implicitly take into 
account renormalization of the temperature in the sys- 
tem. This comes from the fact that creating vortex- 
antivortex pairs removes the energy from the phonon de- 
grees of freedom. We are going to investigate this issue 
in more detail in a separate publication. 
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